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Abstract. In weak gravitational lensing, the image distortion caused by shear measures the projected tidal grav- 
itational field of the deflecting mass distribution. To lowest order, the shear is proportional to the mean image 
ellipticity. If the image sizes are not small compared to the scale over which the shear varies, higher-order distor- 
tions occur, called flexion. 

For ordinary weak lensing, the observable quantity is not the shear, but the reduced shear, owing to the mass- 
sheet degeneracy. Likewise, the fiexion itself is unobservable. Rather, higher-order image distortions measure the 
reduced flexion, i.e., derivatives of the reduced shear. We derive the corresponding lens equation in terms of the 
reduced flexion and calculate the resulting relation between brightness moments of source and image. Assuming 
an isotropic distribution of source orientations, estimates for the reduced shear and flexion are obtained; these are 
then tested with simulations. In particular, the presence of flexion affects the determination of the reduced shear. 
The results of these simulations yield the amount of bias of the estimators, as a function of the shear and flexion. 
We point out and quantify a fundamental limitation of the flexion formalism, in terms of the product of reduced 
flexion and source size. If this product increases above the derived threshold, multiple images of the source are 
formed locally, and the formalism breaks down. Finally, we show how a general (reduced) flexion field can be 
decomposed into its four components: two of them are due to a shear field, carrying an E- and B-mode in general. 
The other two components do not correspond to a shear field; they can also be split up into corresponding E- and 
B-modes. 

Key words, cosmology - gravitational lensing - large-scale structure of the Universe - galaxies: evolution - galaxies: 
statistics 
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>- 1. Introduction 



Weak gravitational lensing provides a powerful tool for studying the mass distribution of clusters of galaxies as well as 
the large scale structure in the Universe (see Mellier 1999; Bartelmann & Schneider 2001; Refregier 2003; Schneider 
2006; Munshi et al. 2006 for reviews on weak lensing). It has led to constraints on cosmological parameters, such as 
those characterizing structure formation and the mass density of the Universe. 

In weak lensing, one employs the fact that the image ellipticity of a distant source is modified by the tidal 
gravitational field of the intervening matter distribution. Based on the assumption that the orientation of distant 
sources is random, the ellipticity of each image yields an unbiased estimate of the line-of-sight integrated tidal field, 
usually called shear in lensing. The shear thus carries information about the properties of the mass distribution. 
Formally, the shear is described in terms of a first-order expansion of the lens equation, i.e., the locally linearized 
lens equation. This yields a valid description of the mapping from the image to the source sphere, as long as the 
images are small compared to the length-scale on which the shear varies. However, this linear approximation breaks 
down for larger sources, or in regions of the lens plane where the shear varies rapidly. The most visible failure of the 
linearized lens equation is the occurrence of giant arcs, which in most cases correspond actually to multiple images 
of a background source; to model them, the full lens equation needs to be studied. However, there is an intermediate 
regime where the linearized lens equation breaks down, although (locally) no multiple images are formed - the arclets 
regime. Arclets are fairly strongly distorted images of background sources (Fort et al. 1988; Fort & Mellier 1994), 
though they do not correspond to multiple images. 
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Arclets are the most natural application for flexion. Flexion has been introduced by Goldberg & Bacon (2005) 
and Bacon ct al. (2006), and describes the lowest-order deviation of the lens mapping from its linear expansion. It 
corresponds to the derivative of the shear; in combination with a strong shear, this can deform round images into 
arclets, giving rise to images which resemble the shape of a banana. In their original paper, Goldberg & Bacon 
(2005) considered only a single component of flexion which, however, only provides an incomplete description of shear 
derivatives. In Bacon et al. (2006), the need for a second flexion component was recognized. 

In the first part of this paper, we present the general theory of flexion; in contrast to earlier work, we explicitly 
consider the quantities that can be actually observed, by accounting for the mass-sheet degeneracy (Falco et al. 1985; 
Gorenstein et al. 1988). That is, a change of the surface mass density k of the form k ^ Ak + (1 — A) leaves the shape 
of all observed images invariant. In usual weak lensing, this is accounted for by recognizing that not the shear 7 can be 
obtained from observations, but only the reduced shear g = 7/(1 — k) (Schneider & Seitz 1995). The difference of shear 
and reduced shear is typically small, in particular in applications of cosmic shear, since along most lines-of-sight, the 
value of K is very much smaller than unity. In applications of flexion, however, we expect that the surface mass density 
no longer is very small; for instance, arclets occur in the inner parts of clusters where k <; 0.1. Therefore, the difference 
between shear and reduced shear can no longer be neglected. Gradients of the shear are not directly observable; only 
derivatives of the reduced shear are, and thus we define the (reduced) flexion in terms of derivatives of g. In Sect. 12.1] 
we briefly recall the irreducible tensor components which are defined in term of their behavior under rotations of the 
coordinate system. It turns out that a complex notation for these tensor components is very useful. In Sect. 12. 2] we 
expand the lens equation to second order, before deriving the corresponding lens equation (and relation for the local 
Jacobian) which is invariant under mass-sheet transformations. The second-order term in this lens equation is fully 
described by our reduced flexion components Gi and G3. 

As is known from usual weak lensing studies, a measured shear is not necessarily accounted for by an (equivalent) 
surface mass density. Since the shear is a two-component quantity, it has one degree of freedom more than the k field. 
Therefore, shear fields are decomposed into E- and B-modes (Crittenden et al. 2002; Schneider et al. 2002), where 
the former are due to a k field, whereas the latter describes the remaining ("curl") part. A similar situation occurs in 
flexion, which has four components. Therefore, in Sect. [3] we consider the decomposition of a general flexion field into 
contributions due to the gradient of the shear and those not related to the shear field. The former one can then be 
further subdivided into flexion resulting from an E- and B-mode shear field. We carry out this decomposition for the 
flexion as well as for the reduced flexion. 

In Sect. |4] we then deflne brightness moments of sources and images and derive the transformation laws between 
them. This approach is very similar to the HOLICs approach developed by Okura et al. (2006) and later considered 
by Goldberg & Leonard (2007), except that we explicitly write all relations in terms of the reduced shear and the 
reduced flexion. Generalizing the usual assumption that the expectation value of the source ellipticity is zero - due 
to the phase averaging over source orientations - to the expectation values of all source shape parameters which are 
not invariant under coordinate rotations (as appropriate for a statistically isotropic Universe), we obtain in Sect. [5] 
estimates for the reduced shear and reduced flexion in terms of the brightness moments of the images. In Sect. [6] we 
perform a number of numerical experiments to test the validity of our approach and the accuracy of the estimators 
derived. In particular, we point out that there is a fundamental limit where the theory of flexion has to break down 
- the second-order lens equation is non-linear and will in general have critical curves, leading to multiple images of 
the source (or parts of it). If the source is cut by a caustic, different parts of it will have different numbers of images, 
and the assumption of random source orientation (which underlies all weak lens applications) will break down - the 
caustic introduces a preferred orientation into the source plane. In appendix B we provide a full classification of the 
critical curves of the second-order lens equation and use these results in order to obtain the maximum source size (for 
given values of the reduced flexion) for which the flexion concept still makes sense. We discuss our results in Sect.[71 

2. Complex lensing notation 

Like in many other instances in weak lensing, flexion is best described by using complex notation, which we shall 
briefly introduce next and which will be used for vectors and tensor components throughout this paper. 

2.1. Irreducible tensor components 

For a two-dimensional vector a; = {xi,X2), we define the complex number a; = Xi+ix2. Under rotations of the coordinate 
system by an angle ip, x gets multiplied by the phase factor e~"^. For a tensor of second rank, whose Cartesian 
components are Qij, we define the complex numbers Q2 = Qii — Q22 +i((9i2 + Q21) and Qo — Qii + Q22 + i(Qi2 — Q2i)- 
A rotation of the coordinate systems by an angle ip multiplies Q2 by the phase factor e^^"^, whereas Qo remains 
unchanged. This is most easily seen by considering that the prototype of a second rank tensor is Qij = XiHj, where 
X and y are vectors; the foregoing statements are then obtained by noting that the complex numbers xy and x*y 
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are multiplied by e~'^^^ and 1, respectively, under coordinate rotations. According to this transformation behavior, we 
shall loosely speak about Qq as a spin-0 quantity, whereas x and Q2 are spin-1 and spin-2 quantities, respectively. 
We shall be dealing only with totally symmetric tensors. If Qij is symmetric, then 

Q2 ■■= Qii - Q22 + 2iQi2 ; Qo := Qii + Q22 ■ (1) 
If TyTj is a symmetric third-rank tensor, we define its spin-3 and spin-1 components as 

T3 := Till — 3Ti22 + i (3rii2 — T'222) ; Ti ;= Tm + ri22 + i (T112 + T222) • (2) 

Furthermore, if Fijki denotes a symmetric fourth-rank tensor, we decompose it into its spin-4, spin-2 and spin-0 
components, respectively, 

F4 := Fun— 6Fii22+-F2222+4i (F1112 — -^1222) ; F2 -Fllll— ^2222+21 (F1112 -f i^l222) ; Fq Fiiii+2i^ii22+-F2222 -(3) 

Apart from notational simplicity, the complex lensing notation provides a check for the validity of equations. In a valid 
equation, each term has to have the same spin. The product of a spin-m and a spin-n quantity has spin m + n. The 
complex conjugate of a spin-n quantity has spin — n. 



2.2. Second-order expansion of the local lens equation 

In weak lensing, the lens equation is linearized locally by writing the relative source coordinate (3 in terms of the 
image position 6 as j3i ^ di ^ ip^ijOj, where is the deflection potential, indices separated by a comma denote partial 
derivatives with respect to 9i, and summation over repeated indices is implied. Note that the form of this equation 
implies that the origin of the lens plane, = 0, is mapped onto the origin of the source plane. The surface mass 
density k and the complex shear 7 at the origin are given in terms of the deflection potential, k = {ip^n + V'.22)/2, 
7 = (V',11 — V", 22)72 -l- iV',12, being spin-0 and spin-2 fields, respectively. In our complex notation, the locally linearized 
lens equation reads 

/3 = (1 - k)6' - 76I* . (4) 

We next generalize this result to a second-order local expansion of the lens equation, which in Cartesian coordinates 
reads (3i = 9i — ip,ijOj — 'ip^ijk0j9k/2. The third-order derivatives of ip are related to the gradient of k and 7. To write 
these derivatives also in complex form, we define the differential operators 

^ d . d d . d 

^''■=w,^'W2'^ ^^'■^wr'W2- 

The differential operator Vc turns a spin-n field into a spin-(r7, -f- 1) field, whereas V* reduces the spin by one unit. 
One finds, for example, 

VcK ^ + V',122 + i (■(/',112 + l/'.222)] ; Vc7 = ^ [l/'ail - 3?/',122 + i(3l/'412 - "0.222)] ; V*7 = VcK , (6) 

and we recognize the combinations of third derivatives of ip which form the spin-1 and spin-3 combinations defined in 
The final relation in ^ is the relation between first derivatives of k and 7 found by Kaiser (1995), here expressed 
in compact form. It expresses the fact that the third-order derivatives of the deflection potential can be summarized 
in the spin-3 field G = Vc7 and the spin-1 field ^ = V*7, where we introduced the usual notation for the two flexion 
quantities. The second-order lens equation in our complex notation then reads 

/? = (1 - k)6I - 76I* - ij^* 9^ - ]^T99* - {9*f . (7) 

Since this is no longer a linear equation, a source at (3 may have more than one image. In fact, up to four images of 
a source can be obtained, as can be seen for the special case of 7 = = and by placing the source at /3 = 0. In 
this case, if we set Q = j^je"^'^, then one solution is = 0, and the other three are 9 = 4(1 — K)/|t/| e"'', with (p — 
(fi — C, + 2tt/3 and ip = ( + Att/S. Of course, the origin for the occurrence of these solutions lies in the fact that ^ is a 
spin-3 quantity. We shall later need the Jacobian determinant det A of this lens equation, which is 

detA = (1 - k)2 - 77* + 6/ • V [(1 - nf - 77*] + 0{e^) 

7*0 + 7-^* 



= (1 - k) - 77* 



(1-^)-^* + ^ V 



(1 - k)J^- 



-0{9^), (8) 



where the flrst expression is just the flrst-order Taylor expansion of the Jacobian around the origin, and in the second 
step we made use of the relation V — (0V* + 6'*Vc)/2. We point out that ([8]) is not the full Jacobian of the lens 
equation ([7|), but only its flrst-order expansion; the full Jacobian contains quadratic terms in 9. We will return to this 
important issue further below. 
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2.3. Accounting for the mass-sheet degeneracy 

The observables of a gravitational lens system are unchanged if the surface mass density n is transformed as k(6) — > 
k'{6) = Xk{6) + (1 — A) (Gorenstein et al. 1988). In the case of weak lensing, the shape of images is unchanged under 
this transformation (Schneider & Seitz 1995). Because of this mass-sheet degeneracy, not the shear is an observable in 
weak lensing, but only the reduced shear g = 7/(1 — n). In fact, since we expect that the most promising applications 
of flexion will come from situations where k is not much smaller than unity, the distinction between shear and reduced 
shear is likely to be more important for flexion than for the usual weak lensing applications. Hence, at best we can 
expect from higher-order shape measurements to obtain an estimate for the reduced shear and its derivatives. For this 
reason, we shall rewrite the foregoing expressions in terms of the reduced shear. 

The mass-sheet transformation is equivalent to an isotropic scaling of the source plane coordinates. Hence, we 
divide ([7]) by (1 — k) to obtain 

/3 = 7-^ ^e~ge*-n0^~ 2^, 69* - VI/3 {e*f with = ^t-^ ; *3 = 777^ • (9) 

We will now express the coefficients in the lens equation (|9]) in terms of the derivatives of the reduced shear, 

_ „ T + gT* ^ _yj G + gT ,^„. 
^i = ^c.9=^^3^; G3 = Ve.g=^^. (10) 

The expression for J- /{\ — k) in terms of the reduced shear and its derivatives has been derived by Kaiser (1995); in 
our notation it reads 

^ -Vcln(l -n)^ ^ ^1 ^ f/^ ' . (11) 



[l-n) ^ ' ' \-gg* 4(1 -.g.g*) 

The expression for the derivative of 7 in terms of the reduced shear can be easily obtained from differentiating the 
definition 7 = (1 — K)g, 

Vc7 _ Q ^ ^^^G^^ 9iG,-gGl) 



(1-k) ^(1-At) 1-55* 4 4(1-55* 

The derivatives Gi,3 of the reduced shear are those quantities we can hope to observe; to distinguish them from T 
and CJ, one might call Gi_3 the reduced flexion. 

The Jacobian determinant det A of the mapping between the image position 6 and the rescaled source position (3 
then becomes 

" (1 - k)2 = ^ ~ f 5 -T] e-r]9 , where 77 = V^5 + — - — = Gi — + — ^ (13) 

is a spin-1 quantity. Again, (|13p is valid only to linear order in 6. Note that a similar equation for the determinant was 
obtained in Okura et al. (2007; their eq. Al), but they consider only the case of \g\ ^ 1; this has also consequences 
for the relations between source and image brightness moments, to be derived further below. 



3. Compatibility relations 

Flexion has a total of four components, namely the real and imaginary parts of J- and Q. A measurement of flexion 
will thus yield four components, and we might ask whether these components are independent. We recall a similar 
situation in shear measurements. The shear has two components; on the other hand, the shear is defined as second 
partial derivatives of the deflection potential, which is a single scalar field. Therefore, the two shear components 
cannot be mutually independent if they are due to a gravitational lensing signal. Of course, the measured shear is 
not guaranteed to satisfy the condition that the two shear components can be derived from a single scalar deflection 
potential, since observational noise or intrinsic alignments of galaxies may affect the measured shear field. Therefore, 
one has introduced the notion of E- and B-modes in shear measurements (Crittenden et al. 2002). The E-mode shear 
is the one that can be written in terms of a deflection potential, whereas the B-mode shear cannot. 

Formally, the E- and B-mode decomposition can be written in terms of a complex deflection potential tp{(^) — 
ip^{0) + iV'^(^) and a complex surface mass density k = + ik^ (Schneider et al. 2002). Each component of ip 
satisfies its own Poisson equation, V^tA^ — 2k^, V^V^ — 2kP. Making use of this decomposition, the shear becomes 



7 = 71 + 172 = (V',11 - "^.22) /2 + iV',: 



^\2 



1 



(14) 
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The distinction between E- and B-mode shear can be obtained by considering second partial derivatives of the shear 
components. Taking the derivative of (HU, one obtains 

T = = (1/2) (^^'ni + - tii2 - i^%2) + (i/2) + + V'I'iii + ^^22) = " 4 + ^ (4 + ' (15) 

which can be expressed in more compact form as 

= Vc (k^ + iK^) = VcK . (f6) 

A further derivative yields for the components 

T — . T7 _ ,.E ^B . T7 _ K I ,^B . T7 _ , ^B /-i 7\ 

-/"la — ^^,11 — K_i2 , -^1,2 — '^a2 ^ '^,22 I -^2,1 — K,12 + 7 -^2,2 — K,22 + '^,12 • U'7 

However, it is easier to consider directly the complex derivative of JF, from which we obtain 

= v: v:7 = (^^ + i^^) = + T2.2 + 1 (^2,1 - ^1,2) . (is) 

Thus, if the shear field is a pure E-mode field, V*V*7 is real. An imaginary part of V*V*7 is due to a B-mode field. 
This then yields the local distinction between E- and B-mode shear. 

Since the flexion has four components, whereas the lens can be described by a single scalar field, we expect that 
there are three constraint relations a flexion field has to satisfy if it is due to a lensing potential. In fact, even if we 
leave the shear field arbitrary (that is, even if we allow it to be composed of E- and B-modes), then we expect two 
constraint equations, since the flexion fleld has two components more than the shear field. These constraint equations 
are easy to obtain. First, if the flexion field is due to a shear field, then we have 

VcV:7 = V,*Vc7 ^ H:= Vc.F-V,*g==0, (19) 

where we defined the spin-2 quantity Ti. It may describe contributions to the flexion which are not caused by a 
shear field, such as due to noise, intrinsic source alignments or higher-order terms (such as lens-lens coupling) in the 
propagation equation for light bundles. As a spin-2 field, a non-zero Ti can be decomposed into its E- and B-modesQ 
If 7i = 0, then the spin-3 fiexion Q is completely determined by the spin-1 flexion T up to an additive constant, as can 
be best seen in Fourier space, for which ^9]) yields Q{t) = — i7(^) t — {I j i*)T {pj ■ Second, if the flexion field is solely 
caused by a gravitational lens effect, i.e., by a pure E-mode shear field, then V*J-' is real, i.e., 

J^i "^IT - VcJ^* = . (20) 

Thus, fiexion from a pure E-mode shear field is characterized by the three constraint equations 7i = and T\ = 
where the former is a two-component equation. 

Turning now to the reduced flexion, the compatibility equations can be obtained as follows. First, if the flexion is 
due to a shear fleld, we have 

:= VcGi-V:G3 = 0, (21) 

as follows from the definition (jlOp of the two flexion components in terms of the reduced shear. Again, if this equation 
is satisfied, G3 is completely determined by Gi, up to an additive constant. Second, if the flexion is caused by a pure 
E-mode shear, i.e., if the shear is due to a real surface mass density, then we employ the quantity ln(l — k), which is 
real and invariant under mass-sheet transformations, up to an additive constant. Therefore, K2 = — V*Vc ln(l — k) 
must be real. We find: 

[v;Gi-.g(VeGi)*] {G\g* ^gG^Gl-G^Gl-g'GlGt) 



i^2 = v:(^) =v: 



1 - .95 



^ (Gi - G\g) 



(22) 



1-55* (1-55*)' 

so that a flexion coming from an E-mode shear fleld satisfies K2 = K^. 

We point out that the foregoing relation suggests a natural way to use fiexion for finite-field mass reconstructions 
in weak lensing. Seitz & Schneider (2001) formulated the finite-field mass reconstruction from measured reduced shear 
in terms of a von Neumann boundary value problem for K — — ln(l — k), whose solution determines K up to an 
additive constant. The 'source' for V^if was determined by the reduced shear and its derivatives, and is given by (HH). 
In Seitz & Schneider (2001), the derivatives of the reduced shear were obtained by finite differencing of g. If fiexion is 
measured, one can replace the 'source' for if by a weighted sum of the differentiated reduced shear field and the 
combination [K2 + i^l)/-^ of the flexion fleld, with the weights chosen according to the estimated noise properties of 
both contributions. 



^ Let Il(ff) be any spin-2 field, and denote by and the E- and B-mode components of H. They can be obtained from 
H most easily in Fourier space, namely 

as can be best seen by taking the shear field as a prototypical spin-2 field; here, (5 is the phase of the complex wave number I. 
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4. Brightness moments of source and image 

We consider an image of a source, and denote the brightness distribution of the source by Since surface brightness 

is conserved by lensing, the brightness distribution of the image is I{9) ~ P{(3{9)). Since the scaling of the source 
plane is unobservable, we shall only work in the following in terms of the scaled source plane coordinates, and therefore 
drop the hat on f3, as well as on A. 

We define the origin of the image (or lens) plane as the center-of-light of the image under consideration, i.e. we 
require 

J d^eei{e) = o. (23) 

Let F{P) be a function of the source coordinate; we define the operator Mom[i^(/3)] as 

Mom[F{P)] = J d^P F{(3) = J d^O det^(0) F{f3{e)) I{9) w J d^O (1 - gg* - rfO - rjO*) F{f3{e)) 1(9) , (24) 

where here and in the following, we use the linear approximation for detyl. In particular, setting F — 1, one finds that 
Mom[l] = S'o = Jd^priP) = J d^e {1- gg* -7fe-r]e*) 1(9) = {1- gg*) S ^detAoS , (25) 



since first-order moments of the light distribution in the lens plane vanish, due to our choice ([23|) of the coordinate 
system. Here, S is the flux of the lensed image, so that S = Sq/ det^o, as usual, where dety^o is the Jacobian at the 
origin 9 = 0. 

4.1. Centroid shift 

The origin of the coordinates in the source plane is the image of the origin in the lens plane as mapped with the lens 
equation. In particular, this does not coincide with the center-of-light of the source, which is given hy (3 = Mom[/5]/S'o, 
or 

fd'f3 f3P{P) = ^ . / d'9 (1 - gg* - 7j*9 - ^9*) [9 - g9* - *J 0^ _ 2vl/, 99* - {9*)^] I{9) . (26) 
So J S{l-gg*)J 



Expanding the integrand, we note that terms linear in 9 vanish, due to (|23p . Defining the second-order brightness 
moments of the image in the form 

Q2 = \ j A^O 02 1{9) ■ Q^=^ j d^9 9 9* I{9) , (27) 
we obtain for the source centroid shift 

-= 3Gig*-5Gl-2gGl , AgG*, + g^G*, - G^g* - Gi(3 + gg*) ^ , hgG^ - ig^Gl - {I - ^g*)G^ 
4(1-35*) 2(1-35*) 4(1-53*) 

We now write these equations in a more compact form; for this, we define the matrix G by G"^ = (6*3, G|, Gi, G3), 
where the 'T' denotes the transpose of the matrix. Then, 

/3 = BG , (29) 

where the coefficients of B = (5i, 62, ^3, 64) are given by 

, ^ .g^Qo - 9Q2 , ^ 85Q0 - 5Q2 - 352Q2 
' 2(1 -5.g*) ' ' 4(1-55*) 

, 35*Q2-2(3 + 55*)Qo + 55Q^ (355* - l)Q*^ - 25*go 

°3 = 77^ ; 04 = -rr-. ;t • (3U) 

4(1-55) 4(1-55*) 

The centroid shift in the source plane is thus given by the product of the derivatives of the reduced shear (expressed 
by Gi and G3) and the area of the image, which is proportional to Qo and Q2- Of course, since the reduced shear 
and its derivatives are not directly observable, the centroid shift in unobservable as well. To get an order-of-magnitude 
estimate of /?, we assume that the source has a linear angular size Og, consider the reduced shear to be of order unity, 
and let 0c be the angular scale on which the reduced shear varies. Then, 

Gn = ; Q„ = O {Ql) ^ /3 = O (^^^ . (31) 
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4.2. Transformation of second-order brightness moments 

Next we consider the second-order brightness moments of the source, defined as Q\ = Mom[(/3 — /3)^]/S'o = 
Mom[/32]/S'o - ^2 and Ql = Moni[(/3 - /3)(/3 - ^)*]/So = Mom[/3/3*]/S'o - f3f3*. By defining the third-order brightness 
moments of the image through 

T3 = ^ J d^O 9^ I{9) ; Ti = ^ J d^O 6^ 9* I{9) , (32) 



we obtain 

Ql = Q2- 2gQo + g^Ql + ""V, """ Tg + """^ ^' ' J "'^ ' "' ^i 



2(1-55) 2(1-55*) 



4(1 - .95*) 

, 2g*^G3 + (11 + 3gff*).9*Gi - (7 + 9ffg*)Gt - (1 + 3.9ff*)gGg 

4(1 -gg*) 

_^ 2g^G5 + (11 + 3gg*)gGl - (1 + 3gg*)g*G3 - (7 + 9gg*)Gi 6gGi - 4g^Gt - (1 - 3gg*)G3 ^, ,„ 

+ i(wr^ "^^^ i(T^^r^ "^^"^^ ^''^ 

Note that Qq is real. In a more compact notation, ([55]) reads 

- Q2 - 2gQo + g'Qa + AG - ,9^ (35) 

where the matrix A = (ai, 02, 03, 04) has coefficients 

^ -g^Ti* + 2g^ri - gTg ^ -3g^r3* + g(7 + gg*)T^ - (4 + 3gg*)ri + 2g*r3 
2(1 -gg*) ' 2(1 -gg*) 

2g3r3* - 7g^T* + 8gri - 3r3 g(l - 2gg*)T; - (1 - 3gg*)T* - g*T, 

"2 = ;t71 ; '^4 = :^7^ ;t ■ (,3Dj 

2(1 -gg*) 2(1 -gg*) 

4.3. Transformation of third-order brightness moments 

We now define the third-order brightness moments of the source, separated into a spin-3 and a spin-1 component, 
^ Mom[(^ ^Pf]^ Mom[g3] _ ^ _ Mom^l ^ 3^. _ ^3 , _ s^g^ _ ^3 ^ (37) 

OO "JO "JO -JO "JO 

where we used that Mom[l3'^]/ So = Ql+f3^ and Mom[P(3*]/ Sq = Qo + Pf^*- Similarly, we obtain 
Defining the fourth-order brightness moments of the image by 

Po^^J d^O {99*fl{9) ■ F2 = ^Jd^9 9^9* I{9) ; Fi = ^ j 6^9 9^ I{9) , (39) 

where Fn is a spin-n quantity, we obtain for the third-order moments of the source: 

= T + C G + 0(^3) , (40) 
where the matrix is defined by its transpose T^-'^ = (TJ**, *, Tf , TJ*). The elements of r are 

Ti = T* - 3g*T* + 3g*^Ti ~ g*^n ; T2 = -gr3* + (1 + 2gg*)r* - g*(2 + gg*)ri + g*2T3 ; = r* ; = t* , (41) 

where the last two relations are obvious. The 4x4 matrix C is given explicitly in Appendix|3 each of its elements 
consists of a sum of terms proportional to fourth-order brightness moments, Fn, and terms proportional to squares 
of second-order brightness moments. Okura et al. (2007) and Goldberg & Leonard (2007) have derived expressions 
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similar to (j40p . though using a number of simplifying assumptions (such as j^l <C 1) and (in the latter paper), not 
considering the reduced flexion. 

We will now consider the order-of-magnitudes of the various terms appearing in (|35p and (j40p . Assuming that the 
third-order moments of the sources are small, then the third-order moments of the image are given by the product of 
C and G. With G = O (l/Gc) and C = O (G^), we find that T = (e^/Gc) = O (G^) {ejQ^). To get an estimate of 
the size of the various terms in (j35p . we note that the first three terms on the right-hand side (those proportional to 
the Q„) are of order O (O^), whereas AG = O (ef/Oc) O (l/Qc) and jS^ = O (ef/O^). Hence, the last two terms are 
of equal magnitude in general, each of them being smaller than the first three terms by a factor (0s/6c)^- Only if the 
source is of the same order as the scale over which the reduced shear varies do the last two terms in ([551) contribute. 
In ((40)) . we have neglected the terms f3^, since they are two powers of (0s/Oc) smaller than the terms written down. 



5. Shear and flexion estimates 

5.1. Estimate of the reduced shear 

We see that (HUl) is a linear equation for G, which can thus be solved, 

G = (T" - r) . (42) 
Inserting this into ([55)1 then yields 

Ql = Q2- 2gQo + g^Ql -f A (T^ - r) - (BG)' . (43) 

We are thus left with a single complex equation for 5, which contains the observable brightness moments of the image, 
as well as the unobservable brightness moments of the source. This equation can be used to estimate the reduced shear 
if we make assumptions concerning the properties of the source brightness moments. We assume that the sources are 
oriented randomly, which implies that all quantities with spin unequal zero have a vanishing expectation value. That 
is, we set Q\ = 0, T*' = 0, to arrive at 

Q2-2gQo + g'Q;-AC-V+(BC-V)'-:r(.g) , (44) 

where we have indicated that the right-hand side depends on the reduced shear (in fact it does so in a very complex 
manner). However, since we have argued above that the terms on the left-hand are much larger than those on the 
right-hand side, an iterative solution of this equation is suggested. Assume the right-hand side is given, then we get 
the solutions 



is the complex ellipticity of the image. Obviously, there are two solutions g for a given value of Y . This situation is 
similar to that of 'ordinary' weak lensing, where this ambiguity also occurs: as shown by Schneider and Seitz (1995), 
from shape measurements of background galaxies, and cannot distinguish locally between an estimate g and l/g* = 
gl\g^ ■ The same occurs here; we therefore assume that we pick one of the two solutions, say the one corresponding 
to the '— ' sign; this then yields for small shear g « x/2. It should be stressed that flexion impacts the determination 
of shear from the second-order brightness moments, due to its impact on higher-order brightness moments; hence, in 
general the determination of shear and flexion are coupled. 

We start the iteration by setting Yq = 0. This yields a first-order solution for the estimate of g, 



X 




(1 - v/T~W^) . (46) 

We then use the iteration equations 

(47) 

This procedure converges quickly to one of the two solutions (5, Gi, 6*3); the other solution is obtained by taking the 
sign in the above equations. 

Of course, our approach of setting Ql — yields a biased estimator for g; this is true even in the absence of flexion 
(e.g., Schneider & Seitz 1995). The reason is that, although the expectation value of Ql vanishes, the resulting estimator 
for g is a non-linear function of =^ Q2/Q0 ^^'^ thus biased. The bias depends on the ellipticity distribution of the 
sources. It should be stressed, however, that a modified definition of image ellipticity exist such that its expectation 
value is an unbiased estimate of the reduced shear (Seitz & Schneider 1997). 
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5.2. Estimates for the reduced flexion 

The flexion estimator is given by (|42p . Since the matrix C contains many terms, this is a fairly complicated equation 
in general. A simpler estimate is obtained if we assume that the reduced shear is small, \g\ <C 1, in which case the 
matrix C simplifies considerably - see Appendix. Furthermore, if we assume that the brightness moments of spin ^ 
are much smaller than the corresponding ones with spin 0, then we find the simple relations 

Tr « - 5^^i^Gr ; r| « T3 - . (48) 

If we then set the T^^ — 0, as would be true for the expectation value, then we obtain as estimates for the reduced 
flexion 

Gi « ^Ti ; G3 w — T3. (49) 

Thus, the flexion is then given by the third-order brightness moments of the image, divided by a quantity that just 
depends on the size of the image. Similar relations to (|49)) have been given in Goldberg & Leonard (2007), whereas 
Okura et al. (2007) obtain a different expression for Gi. We will check the accuracy of (|49|) in Sect. [6] below. 

A more accurate estimate is obtained if we consider the reduced shear as well as the ratios of non-zero spin 
brightness moments to zero spin moments (such as IQ2/Q0I or 1^2,4/^0 1) to be of order S, and then expand the flexion 
to first order in the (small) parameter S to obtain 

^ _ m ^ 4 [2F* + 3Fog* - 2Qo(2g*Qo + Q*)] ^ 4 [SFpg - 8F2 - AQojgQo - AQ2)] , 4^^4T3* 



9Fo - 1202 9Fo(4g2 _ sPo) 9(3i^o - 40^)2 1 9Fo{4Q^„ ~ 3Fo) ' 

_4r3 8(5f2 - 9QoQ2)ri 

' 3Fo ^ 9Fo(4g2 - 3Fo) ^ 9i^o(40§ - 3Fo) ' ^ ' 



6. Numerical tests of flexion estimators 

In this section we describe some simulations that we have performed in order to test the behavior of the estimators 
given in the previous section. 



6.1. Description of the simulations 

We model the sources as elliptical Gaussians, truncated at three times the scale 'radius' 9s chosen such that the 
area of a source is independent of its ellipticity. The ellipticity of the sources follows a Gaussian distribution, with 
a dispersion of of i? = 0.4 (i.e., we use the same ellipticity distribution as in Schneider & Seitz 1995). However, 
for reasons explained in the next subsection, we truncate the intrinsic ellipticity distribution at \x^\ < 0.9. For each 
source, we map a grid of pixels from the lens plane to the source plane using the lens equation to obtain the brightness 
distribution in the lens plane. From this distribution, the brightness moments of the image are measured. A shift in the 
lens plane coordinates is applied as to satisfy ([23|) . We then apply the shear and flexion estimators described above to 
the resulting brightness moments Qm Tn and Fn- The shear and flexion estimates are then averaged over the Gaussian 
ellipticity distribution of the sources, in particular over their random orientation. 

It should be noted that flexion is a dimensional quantity cx Q~^. As can be checked explicitly from Sect.lH the way 
flexion appears in the equations is always with one order higher in the source (or image) size than the other terms in 
the equations. As an example, we consider (|40p : the left-hand side and the first term on the right-hand side are cx Q^, 
whereas the coefficients of the matrix C cx 0^. This then implies that the accuracy of the flexion estimates does not 
depend on the magnitude of the flexion and the source size individually, but only on the product G„0s. Therefore, 
the following results are quoted always in terms of this product. 



6.2. Multiple images, and the breakdown of flexion 

As we mentioned before, the lens equation ([7|) can give rise to multiple images. As can be seen from the example given 
after ([7]), if the flexion is sufficiently small, all but one of these multiple images will be located at a large distance 
from the origin, and the central image of an extended source will be isolated. In this case, this central, or primary, 
image (the shape of which we measure here) is not crossed by a critical curve, and thus the source is not crossed by 
a caustic. The multiple images at large distances from the origin then result from the low-order Taylor expansion of 
the lens equation, which most likely breaks down at these image positions anyway; hence, these additional images are 
of no relevance. If, however, the flexion becomes sufficiently large - or if the source is large enough - this is no longer 
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Fig. 1. Constrains on the combination of source size and reduced flexion for the validity of the concept of flexion. 
Each curve shows the dividing line between a circular source of limiting isophote Q being cut by a caustic (above the 
curve) or not (below the curve); in the former case, the assumptions underlying the flexion concept break down. The 
different curves in each panel are for different values of g, chosen as g = 0.4, 0.2, 0.1, 0.05, 0, as indicated by different 
line types. Without loss of generality, we choose g to be real and non-negative. The four panels differ in the phase of 
the reduced flexion, as indicated. E.g., in the upper left panel, the phases of Gi, G3 are the same as that of g 



the case, and the multiple images of an extended source will merge. If that happens, the whole method of determining 
shear and derivatives thereof from brightness moments will break down. This can be most easily seen by considering 
the caustic curve cutting the source. Different parts of the source will be mapped onto a different number of image 
points in the lens plane, and the caustic curve introduces a direction into the situation. Hence, the assumption of 
an isotropic orientation of sources can no longer be employed. Mathematically, this can be seen from ([M]) : there, the 
transformation between source and image plane no longer is correct if multiple images do occur. More precisely, the 
transformation between source and image coordinates in the calculation of the brightness moments implicitly assumes 
that within the limiting isophote of the primary image, the lens equation is invcrtible. Owing to what was said above, 
the condition that the central image is isolated (so that locally no multiple images occur) can be expressed solely 
by the products G„8s. These products approximately measure the fractional change of the reduced shear across the 
image of a source. 

In our simulations we can check whether a critical curve crosses our central image, just by controlling the sign of the 
Jacobian determinant (the true one, not the linear approximation eq. llSp . If the source size becomes too large, some 
points in the image will have a negative Jacobian. In the Appendix B, we consider the critical curves and caustics of 
the lens equation Q, which allows us to determine the regions in flexion space where no local multiple imaging occurs. 
Some examples of this arc plotted in Fig.[T] Each panel shows the dividing line between parameter pairs (GiO, G3O) 
for a circular source of limiting isophotal radius 0; below the curves, no local multiple images occur, whereas for 
parameter pairs above the lines, the flexion formalism using moments necessarily breaks down. The different lines in 
each panel correspond to different values of g. The occurrence of critical curves also is the reason why we truncated 
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the intrinsic ellipticity distribution of the sources in the simulations, since in the hmit of jx'^l — > 1, keeping the source 
area fixed, there will be orientation angles for which the source will hit a caustic. 

6.3. Estimates of the accuracy 

We now present some results of our numerical simulation regarding the accuracy with which the reduced shear and 
flexion can be obtained with our moment approach. For given input values of g, Gi and Ga, we either measure 
the brightness moments for a single circular source, or average the results over an ellipticity distribution, as described 
above. It should be noted that we have to deal with a 5-dimensional parameter space, namely the 3 complex parameters 
g, Gi and G3, minus one overall phase that can be chosen, e.g., to make g real and positive. Thus, instead of sampling 
the parameter space comprehensively, we only give a few selected results. 

We start by considering a circular source, and determine the effect of flexion on the determination of the reduced 
shear. The left-hand panel of Fig. [2] shows contours of constant fractional deviation Ag/g, in the flexion parameter 
plane. Here it is assumed that the phase of both flexion components is the same as that of g (as would be the case in 
an axially-symmetric lens potential). Errors of order 5% occur already for y^\Gf + G§|8s ~ 0.03, and the fractional 
error increases approximately linearly with the strength of flexion (or with the source size), although it does not 
scale equally with both flexion components. The reason for this effect has been mentioned before - flexion affects the 
transformation between source and image quadrupole moments, as can be seen in p3p . 

In Fig.[3l we show the expectation value of the reduced flexion components, as a function of the input flexion. 
The expectation value has been determined by averaging over an isotropic ensemble of elliptical sources, as described 
before. The left and right panel show the behavior of the expectation value of Gi and G3, respectively, where the 
other flexion component was set to zero. The dashed curve shows the identity, the plus symbols were obtained by 
using the approximate estimator (|49p . whereas the crosses show the expectation values as obtained by employing the 
full expression (|42|1 . where the corresponding value of g was obtained by the iterative process described in Sect.O It 
is reassuring that the expectation value closely traces the input value, i.e., that the estimates have a fairly small bias. 
Furthermore, we see that the approximate estimator (HS)) performs remarkably well. It is seen that the estimates for 
G3 behave better than those for Gi . This can also be seen from the right-hand panel of Fig.[51 where we plot contours 
of constant fractional error 



where AG„ is the deviation of the estimate of G„ from its input value. For simplicity, we have assumed a circular 
source. We see that the accuracy decreases much faster with increasing Gi than with increasing G3. The reason for 
that may be related to the fact that the estimator of Gi is more strongly affected by the non-linearity of the equations, 
as can also be seen in (150]) . 

7. Conclusions and further work 

In this paper, we have studied the effect of flexion in weak gravitational Icnsing. The main results are summarized as 
follows: 

— Owing to the mass-sheet degeneracy, flexion itself cannot be determined, but only reduced flexion. We have therefore 
written the second-order lens equation (which contains the derivatives of the reduced shear, i.e., flexion) as well as 
the relations between the brightness moments of source and image strictly in terms of the reduced shear and the 
reduced flexion. 

— We pointed out that a general flexion field can be decomposed into a pair of components which is due to a shear 
field, i.e., its derivatives, and a pair of components not related to shear. The former pair can be further separated 
into flexion due to an E- and B-mode shear, with only the E-mode flexion expected to arise from gravitational 
lensing. For the second pair of components, no physical interpretation is available; if they arise in measurements, 
they are most likely due to noise or intrinsic shape effects of sources. General relations to separate these components 
are given. 

— We derive the relations between low-order brightness moments of source and image, taking into account that 
the presence of flexion leads to a centroid shift, and it also affects the relation between second-order brightness 
moments - and thus the estimate of the reduced shear. Hence, the presence of flexion has an impact on the shear 
measurements. Starting from these moment equations, we obtain approximate estimates for the reduced shear and 
flexion. 

— We point out a limit where the flexion formalism ceases to be valid, namely when the product of source size and 
flexion is sufficiently large that parts of the source are multiply imaged locally, i.e., where a caustic cuts through 




(51) 
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Fig. 2. Accuracy of the estimates for reduced shear and flexion. The left panel shows contour of constant fractional 
error of 5%, 10% and 15%, on the estimate of the reduced shear g, as a function of GiQg, where we chose g = 0.05 as 
input value, and assumed the phases of Gi, G3 to be the same as that of g. The estimate was obtained by solving the 
iteration equations given in Sect.[5l The right panel shows the fractional error levels at 3, 5, and 10% for the reduced 
flexion, as quantified by ([5T|) . where the estimate was obtained again with the iterative procedure. In both cases, we 
assumed circular sources 
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horizontal and vertical axis show Gi0s,« — 1,3. For both panels, we take g = 0.05, and G3 = (Gi = 0) for the left 
(right) panel. The line indicates the input value, the plus symbols show the simplified reduced fiexion estimate (|49p . 
and the crosses result from the full expression of reduced flexion ([^^ . As can be seen from the left-hand panel, the 
full estimator for the reduced flexion yields a more biased result that the approximate expression ([35]); we have not 
found a reasonable explanation for this behavior 



the source. We have quantified this with numerical simulations, and also have given a complete classification of the 
critical curves of the second-order lens equation employed in fiexion studies. 
— We have performed a number of numerical experiments to study the bias of the reduced shear and flexion estimators. 
However, due to the high dimensionality of parameter space, no comprehensive study has been presented here. We 
also point out that only the product of flexion and source size matters in the accuracy of estimates. 

The possible occurrence of critical curves in highly distorted images may provide a serious obstacle to applications 
of flexion. Perhaps the most promising application of flexion measurements are those in regions where the shear field 
varies on small scales, i.e., close to galaxies (and thus can be used for galaxy-galaxy lensing) or in the inner regions 
of clusters. However, if one finds a strongly distorted image of a background galaxy as in the case of the arclet A5 in 
Abell 370 (Fort et al. 1988), how can one be sure that it is not due to a merged double image of the source? Using 
flexion for studying small-scale structure in mass distributions can therefore be affected by the occurrence of multiple 
imaging. 
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Similar to the situation in shear measurements, the moment approach for flexion as presented here must be 
modified in several ways to be applicable to real data. First, brightness moments must be weighted in order not to be 
dominated by the very noisy outer regions of the image. As is known from shear measurements, such a weighting affects 
the relation between source and image brightness moments. Secondly, one needs to account for the effects of a point- 
spread function. Both of these modifications have been successfully achieved for second-order brightness moments by 
Kaiser et al. (1995; see also Luppino & Kaiser 1997). Goldberg & Leonard (2007) consider these effects in the context 
of flexion. It should be noted, though, that their consideration of the PSF effects is restricted to unweighted moments, 
for which these effects are given by a simple convolution. In the case of weighted brightness moments, however, the 
PSF effects are much more subtle, and we expect that a formalism in analogy to Kaiser et al. (1995) needs to be 
developed - and that this formalism for higher-order brightness moments will be considerably more difficult than for 
the shear case. 

But even disregarding these complications, the present paper only scratches the surface in investigating estimators 
for reduced flexion and their properties. As mentioned before, the second-order lens equation contains five essential 
parameters. The bias of an estimator for reduced shear and fiexion will depend on these parameters, as well as on the 
intrinsic ellipticity (and higher-order moments) distribution of sources. One might ask whether it is possible to find 
an unbiased flexion estimator, such as was possible to construct for the reduced shear. Unfortunately, we have been 
unable to make analytic progress: even for a circular Gaussian source, the brightness moments of the image cannot 
be calculated analytically. Our ray-tracing algorithm with which we conducted our numerical simulations is almost 
certainly sub-optimal; a more advanced method should be developed to reduce the numerical efforts in calculating 
brightness moments. Beside the bias, it would be interesting to calculate the variance of the various estimators, or 
more precisely, their covariance. 

It may turn out that measurements of fiexion, and PSF corrections, are more conveniently done with shapelets, as 
was originally considered by Goldberg & Bacon (2005), Bacon et al. (2006) and Massey et al. (2006). Even if this turns 
out to be the case (see Leonard et al. 2007 for an application of fiexion measurements in the galaxy cluster A 1689), the 
moment approach provides a more intuitive picture of the effects of fiexion. In addition, the weak lensing community 
has profited substantially from the existence of several different methods to measure shear (see Heymans et al. 2006; 
Massey et al. 2007 for the first results of a comprehensive Shear TEsting Programme, in which these various methods 
are studied and compared); therefore, the development of different techniques for measuring ffexion will certainly be 
of interest once the flexion method will be put to extensive use. 
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Appendix A: The matrix C 

In this Appendix, we list the coefficients of the matrix C which occurs in (|40p : 

4(1 - gg*)Cii = ~2gF; + {%gg* - 3)Fo + 6g*(l - 2gg*)F2 + g*\bgg* - 3)^4 

+ QgQlQo ~ y^9g*Ql + (3 - %g*)QlQ2 + 6g*(4gg* - l)QoQ2 + 3<?*2(i _ ■igg*)Q\ 
4(1 - <75*)Ci2 = ^gFX - 2(5 + %gg*)F* + 9g*(3 + .gg*)Fo - 25*^(12 + gg*)F^ + 75*^^4 - ^gQ*^ + 6(3 -I- ^gg*)QlQ^ 

~ 12.g*(3 + gg*)Ql ~ 3g*(5 -h ■igg*)QlQ2 + 63*^(8 + gg*)Q^Q^ - \hg*^Q\ 
4(1 - <75*)C^i3 = -7F4* + mg*F* - 36g*^Fo + 22g*^F2 ~ 5g*^Fi 

+ IbQf - bAg*QlQo + A:Sg*^Ql + 24g*^Q;Q2 - 42g*^QoQ2 + 9.g**Q^ 
4(1 - gg*)Cii = -2g*Fl + &g*^F; - &g*^FQ + 2g*^F2 

+ Qg*Qf - 18g*^Q*2Qo + l2g*^Ql + Q9*^Q*2Q2 - 6g*^QoQ2 
4(1 - gg*)C2i - 2g^F; - 6g^g*Fo + [Agg*{l + gg*) - 2]F2 + 2g*{l - 2gg*)Fi 

- Gg'QlQo + 4g(l + 2gg*)Ql + Qg^g*QlQ2 + [2 - Agg*{?, + 2gg*)]Q^Q2 + 2g*{Agg* - l)Ql 
4(1 - gg*)C22 = -^g^Fl + 2g{l -f Agg*)F* - 3[3 -f gg*[% -f gg*)\F,, + 25*(8 -I- hgg*)F2 - lg*^F^ + %^Qf 

-2.g(13 + 8.g5*) Q^Qo + 4[3 -f gg*{% + gg*)\Ql -f [5 + 35* (16 + 3.g.g*)]Q;Q2 - 2g*(16 -I- \\gg*)Q^Q2 + \WQ\ 
4(1 - gg*)G2z = IgF^ - 2(4 -f Qgg*)F; + 3g*(7 + bgg*)F^ - 25*^(9 + 2gg*)F2 + bg*^F^ - l^gQ*^ 

+ (16 + 385<?*)g^go - 45*(7 + hgg*)Ql - g*(Vi -f \\gg*)QlQ2 + 23*^(17 + ^gg*)Q^Q^ - 9g*^Ql 
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Fig. A.l. The critical curves (left-hand panel) and caustics (right-hand panel) of the lens equation ([9]) for the cases 
of hyperbolic critical curves, as described in Sect. IB. 2l The parameters chosen here are g = 0.05, Gi = 0.07 -I- 0.0151, 
G3 — 0.03 + 0.005i. A circular source is mapped onto two images, as indicated. If the source size were increased, it 
would hit the caustic, the two images would merge, and the flexion concept would break down. The unit of the reduced 
flexion is the inverse of the unit in which coordinates are measured 




Fig. A. 2. Same as Fig. lA.li but for the parabolic case, with parameters g — 0.05, Gi — —0.04, G3 — 0.112 



4(1 - gg*)C2i = (Sgg* - 1)F: - Qgg*^F; + 'ig*\l + gg*)F^ - 2g*'F2 

+ (1 - 7gg*)Qf + 25*(2 + 7gg*)Q;Qo - Ag*\2 + gg*)Ql - 3g*2(i + gg*)QlQ2 + 6g*'^QoQ2 

The other eight elements follow trivially from the foregoing ones, since the second half of the matrix is just the complex 
conjugate one of the first half, i.e., G44 = G*^, G34 = C21 etc., or in general, G^ = 



Appendix B: Critical curves and caustics 

In this Appendix we consider the critical curves of the lens equation ([9]). For this, we need to derive the full Jacobian, 
which can most easily be obtained from considering 9 and 9* as independent variables, and then use d/d9i — d/d9 + 
d/d9\ d 139-2. = i (9/96* - dld9*), which can be inverted to yield 5/36* = V;/2, dld9* = Vc/2. With these relations, 
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Fig. A.3. Same as Fig.|XTl but for the elliptical case, with parameters g = 0.05, Gi = 0.015 + 0.0351, = 0.19 + 0.1051 



one finds that deiA = {dfi / d9){d(}* / dO*) ~ {dp/d9*){dp* /dO) = (V*/3 Vc/3* - Vc/3 V*/3*) /4. Carrying out these 
derivatives, the Jacobian becomes 

det^ = 1 - gg* - if 9 - r]9* + A*9^ + B99* + A{9*f , (B.l) 
with 

A = 4(*? - *J*3) ; B = 4(*i*J - *3*5) ; r/ = 4*i + 2^^-^ + 2g**3 . (B.2) 

Note that A is a spin-2 quantity, whereas is a real scalar, i.e., has spin-0. In the generic case, the critical curves 
(det^ — 0) are conical sections, which may be degenerate, though. We will now perform a complete classification of 
cases that can occur, as well as to derive the critical curve(s) in parametric form. As we shall see, the type of conical 
section is determined, amongst other parameters, by the discriminant 

A = - 4:AA* . (B.3) 



B.l. Zero discriminant 

We start with the case that A = 0, which implies B'^ = AAA*, or B = ±2\A\. The case A = = B either imphes 
that ^1=0 = ^'3, in which case also 77 = so that no critical curves occur, or that '^3 = ^'^/^'*, for which 77 ^ 
in general. In this latter case, the critical curve is a straight line, satisfying r]*9 + ri9* = 1 — gg* . As can be seen by 
inspection, it reads 

9 = ^--^ + iAr? , -00 < A < 00 . (B.4) 

If ^ ^ 0, the phase of A is defined. Since it is a spin-2 quantity, we write A = \A\ e^"^-*. Furthermore, we introduce 
the rotation 9 = xe"'''*. Then the equation for the critical curve reads 

(x ± x*)^ = v*x + vx* + — , with V = — , (B.5) 

1^1 |A| 

and the sign on the left-hand side of the equation depends on the sign of B, where we used B — ±2|yl.|. The parametric 
form of the critical curve, which takes the form of a parabola, can then be written as 

(i^*-!^) \ 2\A\ ; ' i,.*+iy)\ 2\A\ ; ' 

where the first (second) equation applies for B > {B < 0). Note that the parabola degenerates into a straight line if 
u is real (for i? > 0) or purely imaginary (for B < 0). 
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B.2. Non-zero discriminant 

If A 7^ 0, we can perform a translation to eliminate the linear term in det A. Hence we define 9 = 9o + and choose 
^0 such that terms linear in vanish. We then obtain for 9o and for the critical curve condition 



Br, - 
^0 — 1 



with 



A*!^'^ + BM* + A{^*)^ = C , (B.7) 



C - B,,*-Ai,r-AW + - 1 . ^gA* + ,M + Bf -1 , (B.8) 

where the second step was obtained by inserting the expression for rj in terms of the ^''s, and in the final one we 
defined V as the expression in the parenthesis. 

As the first case, we consider ^ = and B ^0 (the case A — Q — B was treated above), which implies that = 
and B = —A'^-^^'^ < 0. The equation for the critical curve then reduces to -B|i9p = C. Furthermore, A = B^, and 
C = -1. Thus, the critical curve is a circle of radius 1/(2|^'3|) and center 6*0, or 9 = 9o + c'^/(2|^'3|), < A < 27r. 

Wc now consider the case A 0; then the phase (pA of A is defined, as used before. Introducing a rotation by 
defining d — x e"''^ , the equation for the critical curve becomes 

\A\ + (x*)^] + Bxx* = {B + 2\A\) xf + {B - 2\A\) xl = C . (B.9) 

The presence and topology of critical curves now depends on the signs of A and C . We first consider the case C = 0; 
then, if A > 0, no critical curves occur, except for the isolated point a; = 0. If A < 0, the critical curves are two straight 
lines, as can be obtained from (jB.7[) : inserting the ansatz 'd — Ae"-, one obtains e^''^''^'''^^ = (— S ± i\/— A)/(2|A|). 
Thus, the critical curves are parametrized as 



= 6lo + Ae'^-^W— ^^^^ ; -c5o < A < c5o . (B.IO) 



For the case of C 7^ 0, the consideration of (jB.9|) yields the result that for A < 0, the critical curves consist of two 
hyperbolae. From (|B.8[) we see that negative A implies C > 0. Also note that A < implies that 2\A\ — B > 0, 
2\A\ + B > 0. The critical curves then read 




-00 < A < 00 . (B.ll) 



For the other case, A > 0, we find from (jB.SP that C < 0. If S ± 2\A\ > 0, we then see from ()B.9|) that no critical 
curves exist. If B ± 2\A\ < 0, which in particular implies B < Q, the critical curve is an ellipse parametrized as 

e^VAy / cosA . sinA \ „ , „ 

9 = 9o + — ^ , +i , ; 0<A<27r. (B.12 

VA \^-2\A\-B ^2\A\-B) ' - 

This concludes the classification of critical curves of the lens equation © . The caustics are obtained by inserting the 
parametrized form of the critical curves into the lens equation. In order to see whether a critical curves cuts through 
the primary image of a circular source of outer isophotal radius 0, we calculate the minimum value /3niin of |/3(A)| along 
the caustics. If /3min > ©, the image is not cut by a critical curve. For an elliptical critical curve, the maximum source 
size allowed is /?min; these values are plotted in Fig.[TJ In the cases where two critical curves exist (e.g., two straight 
lines or hyperbolae), the situation is slightly more complicated. Consider, e.g., the case of two straight critical curves. 
Only those sections of them that are closer to the origin are relevant for this consideration, since if the primary image 
of the source is not cut by these closer sections of critical curves, it will still be an isolated image; the caustics coming 
from the outer sections of the critical curves correspond to multiply imaged source sections of secondary images. 
Accounting for this complication, the maximum sources size have been obtained, as plotted in Fig.[T] 
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